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Abstract 

Improvements in the modeling of radiation in low density shock waves with direct 
simulation Monte Carlo (DSMC) are the subject of this study. A new scheme to determine 
the relaxation collision numbers for excitation of electronic states is proposed. This scheme 
attempts to move the DSMC programs toward a more detailed modeling of the physics 
and more reliance on available rate data. The new method is compared with the current 
modeling technique and both techniques are compared with available experimental data. 
The differences in the results are evaluated. The test case is based on experimental 
measurements from the AVCO-Everett Research Laboratory electric arc-driven shock tube 
of a normal shock wave in air at 10 km/s and .1 Torr. The new method agrees with the 
available data as well as the results from the earlier scheme and is more easily extrapolated 
to different flow conditions. 
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transition probability 
molecule or atom 
degeneracy 

represents an emitted photon 
Boltzmann constant 
excitation rate coefficient 
molecule, atom or electron 
number density 
number 

collision number 
time 

temperature 
distance coordinate 
energy of an excited state 
average electronic energy 
scattering angle 
wavelength 
collision frequency 
cross section 
mean radiative lifetime 

fraction of state participating in transition 


Subscripts 


a 

absorption 

e 

value for electrons 

eq 

equilibrium 

j,k 

value for states j,k 
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M 


pertaining to particles of type M 
reaction 
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Introduction 

Recent interest in planetary aerobraking, aeroassisted orbital transfer and hypersonic 
flight has led to the need to perform analyses in the transitional and high energy flow regimes. 
At altitudes below the range for free molecule calculations but above the applicability of 
continuum calculations, the direct simulation Monte Carlo (DSMC) technique is particularly 
useful. It can, however, be applied more generally. The DSMC method involves the physical 
simulation of gas particles in a flowfield as they undergo collisions, boundary interactions, 
and chemical reactions. The success of this method has been in its ability to model highly 
nonequilibrium flows such as encountered in reentry shock waves. The intent of the current 
research is to improve the capability of DSMC to model high energy nonequilibrium flows 
involving radiation. 

At reentry velocities in the upper atmosphere, radiation from the shock wave can be 
a significant portion of the overall heat transfer to the vehicle. Accurate predictions of 
the nonequilibrium radiation in shocks is required for efficient design of thermal protection 
systems. Although the ability to model the process of radiation with DSMC has been 
demonstrated by Bird 1 , much of the current modeling is very approximate. Verification of 
the modeling techniques and investigations of alternate techniques are required to improve 
the level of confidence in the radiation predictions. The present study addresses this goal 
through investigating the modeling of collisional electronic excitation. 

The thermal radiation due to bound-bound transitions between electronic states domi- 
nates other forms of radiation under the 10 km/s conditions of the current study. One of the 
main difficulties with modeling this radiation is determining collisional excitation rates for 
the various excited states. The DSMC method uses a collision number which represents the 
fraction of collisions between certain species which lead to excitation of electronic states for 
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one or both of the particles. Because of the incompleteness of the data available on excita- 
tion cross sections, the first selection of these numbers was based on a qualitative knowledge 
of the magnitudes of the various cross sections. The values were further refined to match 
available data from the appropriate temperature regime. A new method has been developed 
which bases these numbers on the available excitation rate data. The rate of production of 
a species can be expressed in terms of a reaction rate and this, in turn, can be expressed in 
terms of reaction time, r. This reaction time, when used in conjunction with the collision 
frequency, v, defines a collision number R = tv. Thus, for given excitation rates, it is 
possible to define excitation numbers which are temperature and species dependent. Such 
an approach enhances DSMC by reducing the empiricism in the modeling. In addition, it 
makes it possible to compare with a wider range of alternate calculations and experimental 
data. 

A sample problem of a 10 km/s shock wave in air at .1 Torr has been used to evaluate the 
newly developed method. This test case was chosen because it represents the conditions of a 
shock tube investigation which was performed by Allen, Rose and Camm in the early 1960’s 2 
and repeated recently by Sharma, Gillespie and Meyer 3 . The calculations are performed 
with a one-dimensional standing shock wave DSMC program 4 . Differences between the 
methods are evaluated. Wherever possible, the results are compared with the available 
measurements. 


The DSMC Method 

The direct simulation Monte Carlo (DSMC) method, developed by Professor Graeme 
Bird 5 , uses a direct physical simulation approach to solving the flow of low density gases. 
Because the molecular structure of a gas at low densities must be accounted for in the 
solution process, the problem lends itself to direct physical description. The current 
applications of the DSMC method involve the modeling of internal energy, nonequilibrium 
chemistry, ionizing reactions, and radiation processes. The modeling of thermal radiation 
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is discussed in detail in a later section. In most of the other particulars, the DSMC code 
used in the current predictions is identical to the code which has been well documented by 
Bird 1 . 

Certain modifications which have been developed by the authors are employed in the 
calculations presented. The differences pertain to the modeling of ionization and electric 
field effects. Complete details of these modifications and their effect on the solution can be 
found in reference 6. The version of the DSMC code which was used for all the computational 
results presented in this paper includes the above modifications. Therefore, when results 
are presented as “Bird’s method” or the “new method” these phrases should be interpreted 
to refer to the differences in radiation modeling only. 

Thermal Radiation 

In a shock wave with partial ionization, there is also electronic excitation and accom- 
panying thermal radiation. The radiation from bound-bound transitions between electronic 
states is known to be significant in 10 km/s flows and is the only form of radiation consid- 
ered in the current study. For more highly energetic flows, the method should be modified 
to include bound-free and free-free transitions. Because of the large number of radiative 
states and because a significant amount of the radiation comes from minor species and 
sparsely populated states, the radiation modeling is handled with a phenomenological mod- 
eling technique. Unlike the procedures for the rotational and vibrational energy, in which 
each molecule is assigned a single state, each excited particle is assigned a distribution over 
all the available electronic states. This allows the use of a practical number of simulated 
particles without experiencing the unacceptable statistical scatter which would result if each 
were assigned a specific excited state and each excitation reaction were modeled directly. 
Since each particle in a DSMC simulation represents a very large number of real particles, 
it is not inconsistent with Monte Carlo procedures for a simulated particle to represent 
a population of real particles with a distribution of excited states. Except as relates to 
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the calculation of collisional excitation numbers, the radiation calculations presented herein 
employ the method of reference 1. A brief description follows. 

The sampled electronic energy, e e , can be written as an average over the electronic 
states considered. 

^ = (£e J JV.)/JV, n = J2n, ( 1 ) 

From this it is concluded that a successful event excites the participating simulated particles 
in accordance with a distribution resulting from the above equation. Therefore, 

N . /N = 9 j ex P( — e j ! kT) 

3 Y.9j ex P(- e j/kT) 

The sum is taken over all energies of electronic states below the dissociation or ionization 
energy. The distribution of states for a particular simulated particle is essentially a 
Boltzmann distribution, but it depends on the energy of the specific collision which produced 
the excitation. There is no reason to suppose that the distribution of electronic states in 
the macroscopic sense is Boltzmann. The list of states considered for molecular and atomic 
species is given in reference 1 and explained below. This list represents an attempt to include 
those transitions which are most significant in a 10 km/s shock wave. While some of the 
differences between the experiment and calculations may result from omissions in this set, 
such investigations are beyond the scope of this study. The information needed to complete 
the above procedure is the relaxation number for electronic excitation. Determination of 
the correct relaxation number is the subject of the current investigation. 

The molecular band system is the same as that employed by Park 7 and involves the 
electronic states of molecular oxygen, neutral and ionized molecular nitrogen, and nitric 
oxide. The molecular band transitions considered and the mean time to spontaneous 
emission are listed in Table 1. The actual time to emission in the simulation is assumed to 
be exponentially distributed about this mean time. One restriction in modeling radiation 
is that the collision routine time step must be less than the minimum radiative lifetime 1 . 
As can be seen from Table 1, some of the lifetimes are quite small. Thus, when radiation is 
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included, a small simulation time step may be required and long calculation times result. 


Table 1. Molecular band system 


Band 

States 

r (s) 

A (microns) 

(ev) 

N 2 , 1 + 

3^2 

1.1 x 10“ 5 

1.06 

1.17 

N 2 , 2+ 

5^3 

2.7 x 10 “ 8 

0.34 

3.65 

o 2 ,s-r 

5 -> 1 

8.2 x 10“ 9 

0.20 

6.20 

NO, (3 

3 -> 1 

6.7 x 10 “ 7 

0.22 

5.64 

iVO, 7 

2 -> 1 

1.16 x 10“ 7 

0.23 

5.39 

iV+,1- 

3 -> 1 

6.7 x 10 “ 8 

0.39 

3.18 


The states for atomic radiation were combined to form a manageable number of groups 1 . 
Therefore, the radiative transitions are also grouped and are referred to by number rather 
than spectroscopic code. As a consequence of the grouping, each transition generally involves 
only a fraction of the states in the upper group. This fraction, <f>, is included in the tables 
of radiative transitions. These tables are Tables 2 and 3 for atomic oxygen and nitrogen, 
respectively. 


Table 2. Radiative transitions for 0 


Number 

Transition 


r (s) 

A (microns) 

(ev) 

1 

5^4 

0.6 

3.3 x 10 “ 8 

0.83 

1.49 

2 

6 -> 1 

0.1 

2.5 x 10 “ 8 

0.103 

12.0 

3 

6^4 

0.01 

2 . x 10“ 6 

0.45 

2.76 

4 

6^5 

0.7 

2.5 x 10 “ 8 

0.99 

1.25 

5 

^ — 1 
T 

r- 

0.1 

0.4 x 10 “ 8 

0.099 

12.5 

6 

7^5 

0.4 

1 . x 10“ 7 

0.64 

1.94 

7 

T 

00 

0.1 

5. x 10 “ 7 

0.55 

2.25 
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Table 3. Radiative transitions for N 


Number 

Transition 


r (s) 

A (microns) 

(ev) 

1 

4 -> 1 

0.5 

0.5 x 10“ 8 

0.117 

10.6 

2 

4^2 

0.25 

0.2 x 10“ 8 

0.152 

8.16 

3 

4^3 

0.25 

0.5 x 10“ 8 

0.178 

6.97 

4 

5^2 

0.15 

0.2 x 10“ 8 

0.129 

9.61 

5 

5^4 

0.08 

6. x 10“ 8 

0.907 

1.37 

6 

6^2 

0.1 

00 

1 

o 

■\ — 1 

X 

^ — 1 

0.117 

10.6 

7 

6^3 

0.1 

oo 

1 

o 

^ — 1 

X 

^ — 1 

0.132 

9.39 

8 

6^5 

0.7 

5. x 10“ 8 

1.19 

1.04 

9 

7^4 

0.1 

qO 

1 

o 

^ — 1 

X 

^ — 1 

0.45 

2.76 

10 

8^2 

0.005 

3. x 10“ 8 

0.108 

11.5 

11 

CO 

T 

00 

0.003 

5. x 10“ 8 

0.121 

10.2 

12 

8^4 

0.008 

2.5 x 10“ 7 

0.38 

3.26 

13 

T 

00 

0.04 

1. x 10“ 7 

0.65 

1.91 


A simple physical model has been used in the calculations to estimate the probable 
effects of absorption in the flow. The simple model is thought to be adequate for the visible 
and infrared portions of the spectrum because the gas is essentially transparent to the 
radiation over most of this wavelength range. The transparency assumption looses validity 
in the ultraviolet and the calculations in this wavelength range may be degraded as a result. 
If the number density of absorbing particles is n a in a cell of width Ax, the probability of 
absorption of a photon moving at an angle of inclination 6 to the axis in that cell is 

Axn a cr a / cosd (3) 

In this simple model, a a is assumed to be constant. Each time a radiation event occurs, 
the angle 6 is chosen such that all directions are equally possible and the trajectory of the 
photon is followed until it is absorbed in the flow, hits the surface, or exits from the flow. 
For the present calculation, the energy of the absorbed radiation is not put back into the 
flow. The effect of neglecting this energy will be small because the amount of absorption is 
small. However, absorption in the ultraviolet portion of the spectrum should be the subject 



of further investigation. 


Collisional Excitation Numbers 

There are not much available data on the excitation cross sections of the species involved 
in a real air model. Therefore, the determination of the ratio of collision cross section to 
excitation cross section, which yields the collision number, is not an easy task. Separate 
collision numbers are needed for collisions of each species with neutrals, ions, and electrons. 
A primarily qualitative method was adopted by Bird to demonstrate the capability of DSMC 
to predict nonequilibrium radiation. This method is outlined below. 

The electron-ion and electron-neutral elastic cross sections are of the order of 10 -15 
cm 2 and the data that are available suggest that electron impact excitation cross sections 
are on the order of 10 -16 cm 2 . This indicates a relaxation collision number of about 10. 
Bird argued that, because a successful event was assumed to yield a distribution of states 
instead of a single state, a reduction in the cross sections (or an increase in the collision 
number) was necessary. Table 4 shows the relaxation collision numbers which were used for 
electron and ion impact reactions. The values were increased by a factor of 10 for collisions 
with neutrals. 


Table 4. Collision numbers for excitation, Bird’s method 


Species 

Number 

o 2 

100 

n 2 

100 

0 

500 

N 

1000 

NO 

100 

Nt 

100 


Because the numbers employed in this table are rather approximate, we believe that 
such a procedure detracts from using DSMC in future radiation calculations. It is also our 
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contention that if we have the correct relaxation numbers for the model then there is no 
need to compensate for replacing a level by a distribution. 

The Proposed Method 

The object of the investigation was to provide a method to calculate the relaxation 
numbers which is based on the existing data and does not introduce empiricism into the 
procedure. This can be accomplished by using the following procedure which depends on 
the excitation rate data. 

The reactions which are responsible for electronic transitions between lower state, j, 
and upper state, k, are 

M + Bj < — >■ M + B k ( 4 ) 

and 

Bk < — > Bj + hv ( 5 ) 

The largest rates in these reactions are for collisions with charged particles. However, for 
a slightly ionized gas, neutral particle collisions are also important. The production rate of 
species Bj with number density rij is 

drij/dt = y^ y K M (k,j)n k n M 
k 

— 22 K M (j, k)njn M + '^A(k,j)n k - 22 A (h k ) n j (6) 

k k k 

Multiplying by the energy of level j, ey , and summing it is possible to apply the standard 
Landau and Teller 8 theory to express the result in terms of a local relaxation time, r. 

nde e /dt = d/dt (' ^ e j n j) = 

n M ^2 e 3 S K M(k,j)n k -n M 22 W 22 Km (a k )+22 22 A ( k A)nk- 22 w 22 a (a k ) 

j k j k j k j k 

= n(e eq — ej)/r ( 7 ) 
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If the excitation and deexcitation terms on the left hand side of equation 7 are assumed to 
be of the same order of magnitude as the terms representing the departure of the energy 
from equilibrium on the right hand side, 1/r can be estimated as 

1 j t _ E e j n j(K M (j)n-M + A(j)) ^ 

Eto 

when 

Km(j) = ^2 R M(j,k) 

k 

A (j) = '52 A (j> k ) ( 9 ) 

k 

The above equation, while only an estimate for the appropriate relaxation time, should be an 
improvement over the previous method because it does incorporate the available excitation 
rate data and flowfield property variations. 

The fraction of collisions which lead to excitation is equal to l /tv (or the collision 
number for excitation is tv). The needed rate coefficients for charged particle impact are 
available from Park’s NEQAIR radiation code 9 with the exception of those for molecular 
oxygen and the highest excited level of molecular nitrogen. Those rates were taken from 
Slinker and Ali 10 . 

According to Park, 7 the cross sections for the collisions of neutral heavy particles are 
smaller, approximately by the ratio of the two masses, than those of electron collisions. 
In this study the collision numbers for neutral particles have been determined from 
this relationship, although further investigation of the accuracy of this assumption is 
recommended. 

In order to compare with Table 4, sample collision numbers for the test case are 
given in Table 5. These collision numbers represent the values obtained when the average 
temperature and species densities from the significantly radiating portion of the flowfield 
were used to evaluate the rate coefficient and transition probabilities of equations 6-9. 

Table 5. Collision numbers for excitation, new method 
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Species 

Number 

o 2 

34 

n 2 

4 

0 

1 

N 

1 

NO 

13 

N+ 

1 


While much smaller than the collision numbers used by Bird, these values are in the 
expected range for the neutral particles. The numbers for atoms are smaller than expected 
but within the range of possibility. This may indicate a need for some adjustments to 
the procedure. There is a problem with the value for excitation of N^~ which has very 
high excitation rate coefficients compared with the other species. In this case the value 
was generally less than 1. The program requires the collision number to be 1 or greater. 
Therefore the value for N^~ is set to 1, resulting in a model of the excitation of N^~ which 
is collision limited. Recognizing that the entire procedure is approximate in nature and 
that N+ is known to be a strong radiator which is only slightly affected by nonequilibrium 
phenomena 11 , this adjustment is not expected to significantly degrade the results. 

The Larsen Borgnakke approach which is used for the partitioning of electronic energy is 
also the approach used for the partitioning of the rotational and vibrational internal energy. 
This consistency of approach for the partitioning of internal energy is one of the advantages 
of the DSMC method of modeling nonequilibrium flows. The approach depends on the 
correct calculation of relaxation numbers for determining the percentage of collisions which 
result in the redistribution of internal energy. If these relaxation numbers are determined 
correctly, then it will not matter whether the energy is assigned to a specific internal state, 
as with rotation and vibration, or to a distribution of states, as with electronic excitation. 
Since radiation is determined by the electronic energy, then as long as the average, e, is 
computed correctly, the electronic energy of the system will remain the same whether it 
is assigned to one state or a distribution of states. As a result, when a large sample is 
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considered, the results should be independent of the manner in which e is assigned. Thus, 
it is not necessary to alter the relaxation numbers to account for assuming a distribution of 
states, as was done in Bird’s procedure. 


Results 

Both global radiation values and the spectral distribution are investigated. The recent 
measurements by Sharma, et al. agree well with the earlier shock tube result. Therefore 
direct comparison is restricted to the data of reference 2 at this time. For convenience in 
displaying the results, a shock center is defined as the point at which the density is 6 times 
the free stream density. This is defined as x = 0 in the flowfield plots (Figs. 1 and 3). 
Molecules travel from — x to -\-x on average in the flow. 

The global radiation results (Fig. 1) show large differences between the two collision 
excitation number methods. When Bird’s method is used the intensity of radiation near 
the center of the shock is much lower. The present method predicts larger radiation levels 
overall, but particularly in the nonequilibrium region. 

The general radiation pulse shape from the oscillogram trace during the shock tube 
experiment is pictured in Fig. 2. The calculated emission is plotted against distance but, 
since the wave moves 1 cm in one microsecond, the comparison with the trace is easily 
made. The intensity scale of the trace was not calibrated, so absolute comparisons cannot 
be made. However, the peak nonequilibrium intensity appears to be about one order of 
magnitude greater than the relatively level equilibrium intensity immediately behind the 
shock. The rise in equilibrium intensity as one moves away from the shock was attributed in 
the Allen, et al. paper 2 to the attenuation of the shock velocity as it moves down the tube. 
It is not relevant to the computations. The new method yields a ratio of nonequilibrium 
to equilibrium intensity which is more in line with the oscillogram trace than the results of 
Bird’s method. 

The translational temperature and electron temperature for the two methods are shown 
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in Fig. 3. Peak translational and electron temperatures occur ahead of the shock center, as 
defined by the density ratio, except in the case of the new method electron temperature. The 
translational temperature is slightly lower in the present method, but the most noticeable 
difference is the lower electron temperature in the region just ahead of the shock center. The 
larger amount of radiation corresponds to proportionately less energy in the other modes, 
hence the lower temperatures. 

Nonequilibrium Radiation 

The experimental results for the radiation intensity versus wavelength in the nonequi- 
librium (peak radiation) region are compared in Fig. 4 with the results from the two DSMC 
calculations. The radiant energy per unit volume from each of the transitions listed in 
Tables 1, 2, and 3 was spread evenly over the wavelength range defined by the midpoints 
between this transition and the two neighboring transitions. Both methods have reasonable 
qualitative agreement with the experimental data except for the present method overshoot 
of radiation at .39 microns. In general the present method tends to overpredict the radiation 
in the lower wavelength region while Bird’s method tends to underpredict the radiation in 
the longer wavelengths. The overshoot at .39 microns is a result of the large amount of 
radiation from the IV^~, [1—] state. As mentioned before, the electronic excitation rates for 
this ion were very large compared with the rates for the other species and a collision num- 
ber could not be calculated directly. This uncertainty could be blamed for the calculation 
discrepancy. However, the ion is known to be a strong radiator and a small overestima- 
tion of the concentration of N^~ at this point in the flow could equally well account for the 
radiation overshoot. A detailed comparison of the radiation from separate species and tran- 
sitions is given in Figs. 5, 6, and 7. The contribution to the radiation from the molecular 
species is consistently larger with the new method. The largest increase in the molecular 
contribution is the contribution from N ^~ , as was expected from the radiation overshoot at 
.39 microns evident in Fig. 4. In the region below .2 microns the present method shows 
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a significantly larger contribution from N and 0 atom transitions. While this discrepancy 
might be attributed to the smaller collision numbers for atomic species of the new method, 
differences in chemistry or electron energy could also be a factor. Unfortunately, not only 
is the mechanism for this unclear, but no data exist to compare with the computations in 
this region. 

Equilibrium Radiation 

The spectral distribution of radiation from the near equilibrium region (region of 
relatively constant radiation) behind the shock is compared to the shock tube data for the 
two calculations in Fig. 8. Detailed plots of the contributions from the molecular and atomic 
species are plotted in Figs. 9, 10, and 11. The results are similar to the nonequilibrium 
region in that the new method predictions tend to be larger than those of Bird’s method, 
although the magnitude of the disagreement is less in the near equilibrium region. 

Conclusions 

In an effort to improve the modeling of nonequilibrium radiation with DSMC, a new 
method for determining relaxation collision numbers for electronic excitation has been 
proposed. Results for the new model are compared with results obtained using the method 
currently incorporated in Bird’s DSMC programs 1 . These results are also compared with 
data from an electric arc-driven shock tube experiment 2 . 

The determination of the relaxation collision numbers for electronic excitation was 
previously based on a qualitative knowledge of the magnitudes of excitation cross sections 
compared with collision cross sections. The numbers required some adjustment in order 
to obtain agreement with experiment. A new method was proposed which bases the 
determination of these numbers on the available rate data for electronic excitation and values 
of radiative state lifetimes. It has been shown that the relaxation collision numbers obtained 
with this method can be used without adjustment to provide equally acceptable radiation 
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estimates. With the exception of the contribution to the total radiation from the N^~ ion, 
the results obtained with the new method agree as closely with the experimental values as 
those obtained with Bird’s method. The N^~ predictions may be off for a number of reasons. 
The most probable are errors in collision number due to calculation uncertainties and/or 
errors in ion concentration due to uncertainties in the nonequilibrium chemistry. Further 
work is suggested. More detailed modeling of the effect of absorption in the ultraviolet 
portion of the spectrum is also suggested. 

The methods both suffer from the lack of experimental data. These data are required to 
supply some of the modeling parameters. They are also needed in order to evaluate whether 
the methods are adequate for future DSMC applications. However, judging from the data 
that are available, the new modeling technique appears to be a step in the right direction. 
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Captions 


Fig. 1 Global radiation through the shock. 

Fig. 2 Oscillogram records of radiation from shock tube experiment. 

Fig. 3 Temperature through the shock. 

Fig. 4 Nonequilibrium spectral radiation intensity vs. wavelength. 

Fig. 5 Molecular contributions to total nonequilibrium radiation. 

Fig. 6 Contribution of atomic nitrogen to total nonequilibrium radiation. 
Fig. 7 Contribution of atomic oxygen to total nonequilibrium radiation. 
Fig. 8 Equilibrium spectral radiation intensity vs. wavelength. 

Fig. 9 Molecular contributions to total equilibrium radiation. 

Fig. 10 Contribution of atomic nitrogen to total equilibrium radiation. 
Fig. 11 Contribution of atomic oxygen to total equilibrium radiation. 
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Fig. 3 
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Fig. 5 
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